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Abstract. We study the effect of adding to a directed chain of interconnected systems a 
directed feedback from the last element in the chain to the first. The problem is closely related 
to the fundamental question of how a change in network topology may influence the behavior of 
coupled systems. We begin the analysis by investigating a simple linear system. The matrix that 
specifies the system dynamics is the transpose of the network Laplacian matrix, which codes 
the connectivity of the network. Our analysis shows that for any nonzero complex eigenvalue A 
of this matrix, the following inequality holds: |^| < cot A This bound is sharp, as it becomes 
an equality for an eigenvalue of a simple directed cycle with uniform interaction weights. The 
latter has the slowest decay of oscillations among all other network configurations with the same 
number of states. The result is generalized to directed rings and chains of identical nonlinear 
oscillators. For directed rings, a lower bound a c for the connection strengths that guarantees 
asymptotic synchronization is found to follow a similar pattern: a c = i_ cos (27r/n) ■ Numerical 
analysis revealed that, depending on the network size n, multiple dynamic regimes co-exist in 
the state space of the system. In addition to the fully synchronous state a rotating wave solution 
occurs. The effect is observed in networks exceeding a certain critical size. The emergence of a 
rotating wave highlights the importance of long chains and loops in networks of oscillators: the 
larger the size of chains and loops, the more sensitive the network dynamics becomes to removal 
or addition of a single connection. 

Keywords and phrases: coupled systems, reaction network, eigenvalue, synchronization, 
wave solutions 

Mathematics Subject Classification: 34A30, 34D06, 34D45, 92B20, 92B25 


1. Introduction 

A fundamental question in complex networks is how topology influences the overall network behavior. 
This issue is crucial for understanding a range of phenomena in elementary chemical kinetic systems, 
populations of agents, and processes in the neuronal circuits of the human brain [12]. It is well known 
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that a sufficiently stron g d iffusive coupling will lead to globally asymptotically stable synchronization in 
a large class of systems [25[. Some network topologies, moreover, may give rise to partial synchronization 
|lU2u27| , and networks with dynamically changing topologies were shown to exhibit complex multi-stable 
dynamics (see e.g. [j] and references therein). 

Even when network topology is not changing dynamically, its influence on the overall network dy¬ 
namics is well documented. Examples of how the network topology may affect e.g. coherence of network 
dynamics are provided in 0 • The authors showed that shortcuts in otherwise regular lattices significantly 
reduce the critical coupling strength needed for achieving global asymptotic stability of the synchronous 
state. Hence networks with shortcuts can be considered as more efficient than regular ones in term of 
resources spent, such as the total number of connections and their strength, for reaching and maintaining 
synchronous regimes. This aspect of shortcuts appear to be crucial for forming small-world structures 
[14l . ll9l] in evolving networks. Further examples showin g signifi cant dependence of network dynamics on 


the corresponding connectivity graphs can be found in [ 12 l . 2 


Understanding the problem of how the network topology affects its dynamics is a huge theoretical and 
practical challenge if considered in its full generality. Here we will focus on a much simpler objective. In 
particular we will discuss and analyze two basic and extreme topologies which any network will contain 
as a subgraph: a directed chain, and a directed cycle. Not only may these be considered as basic 
building blocks of arbitrary network topologies; many networks can be reduced to chains and cycles as 


well [16l . l28| . Moreover, recent computational studies revealed that cycles could be important on their 


own for sustaining coherent oscillatory network activity [13]. 

We begin our investigation by analysing the dynamics of a system of coupled neutrally stable linear 
equations. The dynamics are essentially governed by a coupling matrix that corresponds to directed 
interconnections in the system. The equations can also be viewed as a model describing the dynamics of 
damped oscillations in kinetic systems, and thus in what follows we refer to it as such. The results are 
provided in Section [2] In Section [3] we consider a more generalized setting, in which the dynamics of each 
individual node is governed by a nonlinear, albeit semi-passive j25j, oscillator. The equations describing 
oscillators in each node are of the FitzHugh-Nagumo (FHN) type 11]. These oscillators are, in turn, 
an adaptation of the van der Pol oscillator [321. We found that for the case of the directed cycle the 
value of critical coupling needed to maintain globally asymptotically stable synchrony is O ^ 1 _ cos ( 2 7 r /ra) ) > 
whereas the synchronization threshold for systems organized into directed chains does not depend on n. 
Moreover, the error dynamics corresponding to the simple directed cycle rapidly becomes underdamped 
for large n; this enables resonances between the dynamics of individual nodes and the coupling dynamics. 
Further numerical analysis show that not only fully synchronous oscillations may occur in these types of 
network, but also a stable rotating wave solution may emerge. These two dynamic regimes co-exist for 
a broad range of coupling strength and for number of systems. Occasionally we observe a shift towards 
prevalence of rotating waves for large enough n. Section [|] contains a discussion of our findings, and 
Section [5] concludes the paper. 


2. Coupled Neutrally Stable Systems 

Consider the following system of linear first-order differential equations: 

P = KP, 

where P = col(pi,p 2 > • ■ ■ ,Pn) € R”, and matrix K = ( kij ) is defined as follows: 

i _ f Qij ■ Qij b if i 7 ^ j ! 

klJ l - E, 


( 2 . 1 ) 


( 2 . 2 ) 


jm, m^i if ^ j- 

Note that K is a Metzler matri/sQ with zero column sums. Off-diagonal elements kij, i ^ j, of the matrix 
K can be viewed as the connection weights between the z-th and the j-th nodes in the network. The 


1 A Metzler matrix is a matrix with non-negative off-diagonal entries 
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matrix K can be related to the Laplacian matrix L (see e.g. i) of an associated directed network in 
which the overall connectivity pattern is the same except for that the direction of all connections is 
altered. The Laplacian for the latter network is thus L = —K T . Note, however, that this relation does 
not necessarily hold for the original network. 

System (12.11) is a commonly used model of first-order kinetics with a finite number of states. In this 
case the variables pi may represent concentration, probability, or population of these states. A more 
detailed discussion and a kinetics interpretation of the model is provided in Section [I] 

Consider the simple simplex 

A n = \P\Pi > 0 = !} • 

An is clearly forward invariant under the dynamics EH) since it preserves non-negativity and obeys the 
“conservation law” JA pi = const. (The latter follows immediately from the fact that K has zero column 
sums.) Thus any solution P(-;to, Pq) of (12.11) starting from P 0 = P(to) € A n remains in A n for all t > to- 
The invariance of A n under eh can be used to prove certain important properties of K and its 
associated system eh- Two examples are presented below. 

— Equilibria. The non-negative vector P* such that KP* = 0 is known as the Perron vector of K, and 
defines an equilibrium of system (EH). The existence of this vector P* can also be deduced from the 
forward invariance of A n . Indeed, as any continuous map : A n —> A n has a fixed point (Brouwer 
fixed point theorem), <P = exp (Kt) has a fixed point in A n for any t > to- If exp (Kt)P* = P* for 
some P* £ A n and sufficiently small t > to, then KP* = 0 because 

exp (Kt)P = P + tKP + o(t 2 ). 

— Eigenvalues of K. It is clear that K has a zero eigenvalue. In fact, Gershgorin’s theorem implies 
that all eigenvalues of K are in the union of closed discs 

A = {AeC|||A-MI < IMI- 


Thus K does not have purely imaginary eigenvalues. This can also be deduced from the forward 
invariance of A n in combination with the assumption of a positive equilibrium P*. We exclude the 
eigenvector corresponding to the zero eigenvalue and consider K on the invariant hyperplane where 
y„ p-i = 0. If K has a purely imaginary eigenvalue A, then there exists a 2D AT-invariant subspace U, 
where K has two conjugated imaginary eigenvalues, A and A = —A. Restriction of exp (Kt) on U is 
a one-parametric group of rotations. For the positive equilibrium P* the intersection ( U + P*) D A n 
is a convex polygon. It is forward invariant with respect to EH because U is invariant, P* is an 
equilibrium and A n is forward invariant. But a polygon on a plane cannot be invariant with respect 
to the one-parametric semigroup of rotations exp(ATi) (t > 0). This contradiction proves the absence 
of purely imaginary eigenvalues. 

The main result of this section is the following theorem: 

Theorem 2.1. For every nonzero eigenvalue A of matrix K 



< 


7T 

cot —. 
n 


(2.3) 


The proof of this theorem can be extracted from the general Dmitriev-Dynkin-Karlelevich theorems 
[9U20I , but the straightforward geometric proof presented below, which makes use of the forward invariance 
of A n under the dynamics EH, seems to be more instructive. 

Proof. Let us assume that system EH has a positive equilibrium P* £ A n (p* > 0 for all* = 1,..., n). 
For this P*, 

J2 qi i p *i yi'l-i'l’'- 

j 3 
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Figure 1. The polygon A is presented as a sequence of vectors x*. The angle pi between 
vectors x. ( and Xj+i and the angles ai and 7 \ of the triangle with sides x,; and x, + i are 
shown. In the Fig., rotation goes anticlockwise, i.e. 3A < 0. In this case, the polygon 
A is invariant with respect to the semigroup exp(t/C) (t > 0 ) if and only if 5 < ai for all 
i = 1 ,..., m, where 5 is the angle between the vector field /Cx and the radius-vector x. 


Systems m without strictly positive equilibria (but with non-negative ones) may be considered as limits 
of those with positive equilibria. 

Let A be a complex eigenvalue of K and let U be a 2D real subspace of the hyperplanc YliPi = 0 that 
corresponds to the pair of complex conjugated eigenvalues, (A, A). Let us select a coordinate system in 
the plane U + P* with the origin at P* such that restriction of K on this plane has the following matrix 


KA -3A 
SA 5RA 


In this coordinate system 


exp (t/C) 


exp(tJiA) cos(tSA) — exp(t!RA) sin(t3A) 
exp(tJiA) sin(t3A) exp(fJiA) cos(t3A) 


The intersection A = (U + P*) fl A n is a polygon. It has no more than n sides because A n has 
n (n — 2)-dimensional faces (each of them is given in A n by an equation pi = 0). For the transversal 
intersections (the generic case) this is obvious. Non-generic situations can be obtained as limits of generic 
cases when the subspace U tends to a non-generic position. This limit of a sequence of polygons cannot 
have more than n sides if the number of sides for every polygon in the sequence does nor exceed n. 

Let the polygon A have m vertices v 7 (m < n). We move the origin to P* and enumerate these vectors 
x, = Vj — P* anticlockwise (Fig [I]). Each pair of vectors x,,x, :+ i (and x m ,xi) form a triangle with the 
angles ai, pi and 7 ,;, where Pi is the angle between x* and x; + i, and p m is the angle between x m and xi. 
The Sine theorem gives -M- = llSiili _M- = -N-. 

Several elementary identities and inequalities hold: 

0 < ai, pi,ji < 77 Pi = 2 tt; a t + Pi + 7 * = n; 

i (2.4) 

^[sinaj = J^sin 7 j (the closeness condition). 


These conditions (12.41) are necessary and sufficient for the existence of a polygon A with these angles 
which is star-shaped with respect to the origin. 

Let us consider anticlockwise rotation (3A < 0, Fig. [TJ). (The case of clockwise rotations differs only 
in notation.) For the angle S between Kyii and x*, sin<5 = —SA, cos S = —3?A and taniJ = 

For each point x € U + P* (x/P*), the straight line {x + e/Cx | e € K} divides the plane U + P* in 
two half-plane (Fig. [Q dotted line). Direct calculation shows that the semi-trajectory (exp(t/C)x 1 1 > 0} 
belongs to the same half-plane as the origin P* does. Therefore, if 6 < ai for all i = 1 ,... ,m then the 
polygon A is forward-invariant with respect to the semigroup exp (i/C) (t > 0). If S > ai for some i then 
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for sufficiently small t > 0 exp (tK.)x.i A because /Cx, is the tangent vector to the semi-trajectory at 
t = 0. Thus, the polygon A is forward-invariant with respect to the semigroup exp(t/C) (t > 0) if and only 
if 8 < Ui for all i = 1 ,..., to. The maximal 5 for which A is still forward-invariant is d max = min,;{a-,;}. 
We have to find the polygon with m < n and the maximal value of minftai}. Let us prove that this is 
a regular polygon with n sides. Let us find the maximizers a^, ft, 7, (i = 1,..., to) for the optimization 
problem: 

min{ay} —> max subject to conditions (12.41) . (2.5) 


The solution of this problem is that all ay are equal. To prove this equality, observe that minftaj} < Th¬ 
under conditions (12.41) (if all a, > \ then the polygonal chain A cannot be closed). Let min^ ay = a. Let 
us substitute in (12.41) the variables cti which take this minimal value by a. The derivative of the left hand 
part of the last condition in (12.41) with respect to a is not zero because a < Assume that there are 
some ay > a. Let us fix the values of ft (i = 1,... ,m). Then 7, is a function of a,, 7y = n — ft — a,. We 
can use the implicit function theorem to increase a by a sufficiently small number £ > 0 and to change 
the non-minimal aj by a small number too, aj 1—>■ aj — 0; 9 = 0(e). Therefore, at the solution of (12.51) all 
a o = a (3 = • • ■ ,m). 

Now, let us prove that for solution of the problem (12.51) all /3, are equal. We exclude 7* from conditions 
(CD} and write /?, + a < n; 0 < /?,, a; 


to log sin a = log sin(/?j + a). 


( 2 . 6 ) 


Let us consider this equality as equation with respect to unknown a. The function log sin x is strictly 
concave on (0,7r). Therefore, for ay e (0,7r) 


1 


\ ' 


log sin — , • 

\ TO. ' 

\ t= 1 


1 


> — > log sin Xi 




and the equality here is possible only if all ay are equal. Let a* £ (0,7t/2) be a solution of (12.61) . If not 
all the values of are equal and we replace /?,; in (12.61) by the average value, (3 = then the value of 
the right hand part of (12.61) increases and sin a* < sin(/3 + a*). If we take all the ft equal then (12.61) 
transforms into elementary trigonometric equation sin a = sin(/3 + a). The solution a of equation (12.61) 
increases when we replace ft by the average value: a > a* because sin a* < sin(/3 + a*), a £ (0,7t/2) 
and sin a monotonically increases on this interval. So, for the maximizers of the conditional optimization 
problem (12.51) all ft = — and ay = 7, = ^ — —. The maximum of a corresponds to the maximum of to. 
Therefore, m = n. Finally, max{(5} = \ — — and 


max 


/Ml = „* 1. 


1 I'KAI / 


□ 


Remark 2.2. It is important to note that the bound given in Theorem 12.II is sharp. Indeed, let K define 
a directed cycle with uniform weights q 7 e.g. 


/-q 0 0--- 0 q\ 

q-q 0 ••• 0 O' 


K = 


0 0 


0 


q-q 0 


V 0 0--- 0 q-q) 

5 


0 


-0 


0 


0 
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The eigenvalues of K are 


A k = -q + qex p 


/ 2nki 

V n 


k = 0 , 1 ,..., n — 1 , 


cf. 


with i = -v/^T the imaginary unit. Thus = cot Note that for large n, 

7 r n 

cot — w —, 
n 7r 

which means that oscillations in a simple cycle with a large number of systems decay very slowly. 

An important consequence of this extremal property of a simple cycle is that not only transients in 
the cycle decay very slowly but also that the overall behavior of transients becomes extremely sensitive 
to perturbations. This, as we show in the next sections, gives rise to resonances and bistabilities if 
neutrally stable nodes in (EH) are replaced with ones exhibiting oscillatory dynamics. As a model of 
nodes with oscillatory activity the classical Fitzhugh-Nagumo 11] system has been chosen. Our choice 
of this system among various alternatives fTU was motivated purely by its simplicity and relevance for 
modelling behavior of neural systems. 


3. Coupled Nonlinear Neural Oscillators 

Consider a network of FitzHugh-Nagumo (FHN) neurons 

(ij= a (yj - fa) 

1 Vj = Vj - 7 Vj ~ Zj + Uj, 

j = 1 , 2 ,...,n with parameters a, /3 ,7 chosen as 

a = B = — , 7=i 
u 100 > v 10 > ' 3 

The FHN neurons interact via diffusive coupling 

n 

Uj = aJ^qjiiyi - yj ) 

1=1 

with constant a G R, a > 0, being the coupling strength. For convenience, let 


(3.1) 


(3.2) 


y= U=[uu n ), x = (y,z), 

and x(-;xo,cr) denote a solution of the coupled system with the coupling strength a and satisfying 
the initial condition x(0) = Xo- The topology of network connections in (13.21) is characterized by the 
adjacency matrix Q with zeros on the main diagonal and entries identical to the values of qij for i ^ j , 
i, j £ {1,..., n}. The matrix Q is now assumed to be a circulant matrix 


Q 


/0 0 ••• 0 1\ 

1 0 ••• 00 

01 : 

: ■■■ '-.0 0 

\0 ••• 0 10 / 


Thus besides assuming the network structure to be a directed ring we have also assumed the interaction 
weights qji to be identical and, without loss of generality, we have set these weights of interaction to 1 . 
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At first glance, the connectivity pattern specified by Q differs from that specified by matrix K in (fA 
Yet, if coupling (13.21) is rewritten in the vector-matrix notation then the following identity holds 

u = cr(Q - I n )y = -crLy. (3.3) 

As remarked before, the network Laplacian matrix 

L = Diag(^ q-ji) - Q = I n - Q 

3*1 

can be related to the matrix K corresponding to the simple cycle in “reverse” direction as L = —K T . 

In what follows we will employ the notions of semi-passivity and strict semi-passivity that have been 
introduced first in 25|. For consistency, we recall these notions below 

Definition 3.1. Consider a system of first-order nonlinear ordinary differential equations 

x = f(t,x,u(t)),y = h(x) 


where f : R x R" xl-> K™ is a continuous and locally Lipschitz function, h : R” —> R is a continuous 
function, and u : R —> R is a continuous function. Let x(•; to, xo, [it]) be a solution of the Cauchy problem 
x(to;to,Xo, [it]) = Xq, and letU C C° be the set of inputs u for which the solution x(-;to,xo, [it]) is defined 
in forward time. 

The system is called semi-passive if there is a non-negative function S : R" —X R+ (a storage function) 
and a function H : R ra —> R such that for each x(-;to,xo, [it]) the following holds for all t > to in the 
domain of this solution definition: 

S(x(t;to,xo, [u])) — S(xq) < / y(r)u(r) — H(x(r, t 0 , x 0 , [it]))dr, 

Jt 0 

where the function H is non-negative outside a ball in R". 

The system is called strictly semi-passive if the function H is strictly positive outside a ball in R™. 


3.1. Boundedness of solutions in the coupled system 

Lemma 3.2. The solutions of the ring network of FHN neurons are ultimately bounded uniformly in Xq, 
er £ R>o- That is, there is a compact set Cl £ R 2 ™ such that for all xq £ R 2ra , er £ R> 0 

lim dist (x(t, xq, er), 12 ) = 0 . 

t—foo 

Proof. We being begin with establishing that the FHN neuron is strictly semi-passive (see also S3). 

Let S(zj,yj ) = \ [a 1 ^ 2 +y 2 ) be the storage function. Then 

S = —H(zj,yj) + yjUj 

with H(zj,yj) = /3zj + y 2 ( 7 y 2 — l). Noticing that 

P z j + Vj (72/j - !) = flz 2 + dy 2 + 7 y A - y 2 - dy 2 = 

^ 2 + d, 2 +(^y 2 -^) 2 -^if (3 ' 4> 

we can conclude that H(zj,yj) is positive for all Zj , y.j such that 

Pz 2 + dy 2 > (3.5) 

7 


© 


© 


0 


0 














© 


’arXive'ms' 


2015/5/22 — 0:42 — page 8 — #8 


© 


© 


u 


© 


A. N. Gorban, N. Jarman, E. Steur, C. van Leeuwen, I. Yu. Tyukin 


Leaders do not look back 


Assigning the value of d in (13.51) as d = /3, ensures that H(zj,yj ) is positive outside the ball 


Za + 


v?< 


(<3+l) 2 

4/?7 


Now consider V(z , y) = S(zi,yi) + . ■ ■ + S(z n , y n ). Then the strict semi-passivity property of the FHN 
neurons implies 

V < -H(z 1 ,yi) - ... - H(z n ,y n ) - cry T Ly. 

Notice that the matrix L + L T is the Laplacian matrix of the undirected ring, which is known to be 
positive semi-definite. Hence 

y T Ly = \y T (L + L T )y > 0, 

and consequently 


V < -H{z l ,yi) - ... - H(z n ,y n ). 

Therefore, setting the value of d in (13.41) equal to a/3 results in 

n 

V < - + Pay] + n^S +111 = -/3aV + 

i=i 

Noticing that the function V is radially unbounded, positive-definite, we invoke the Comparison Lemma 
(see e.g. 0) in order to conclude that solutions of the coupled system are bounded and converge 
asymptotically to a compact set of which the size is independent of the parameter a. □ 

3.2. Sufficient conditions for synchronization 

3.2.1. Directed chain: ”no looking back” 

First we consider the dynamics of two coupled systems in the leader-follower configuration: 

p 1= a( yi -M) (3 . 6) 

\ yi = yi - 72/i - zi 
( z 2 = a(y 2 - Pz 2 ) 

\ V2 = V2 - 72/2 - Z2 + v(yi - Vi)- 


Theorem 3.3. Consider the system of coupled FHN oscillators \3. 61) in which the parameter a is chosen 
so that 

a > 1. 


Then solutions of the system asymptotically synchronize for all values of initial conditions. 


Proof. In accordance with Lemma 15721 solutions of the coupled system exist and are bounded for all t > 0. 
Define 

Z = Z!-Z2, y = y i - 2/2, 


such that 


z = a(y — ffz) 
y = y - i{yl - yl) - z - cry. 

Consider the function 

v = h& 2 + v 2 ), 
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(a) (b) 


Figure 2. Synchronization of two FHN oscillators for a = 1.5. (a) Outputs of FHN 
oscillator 1 (leader, black) and FHN oscillator 2 (follower, red), (b) Synchronization 
output error y := y\ — 1/2 • 


then, using the equality 

(: 2 /i ~ 2/2)(y? - 2 / 1 ) = |( 2 /i - 2/2X ( 3 (j/i + y 2 ) 2 + (2/1 - 2/ 2 ) 2 ) , 

we find 

V = -fiz 2 + (1 - &)y 2 - 4y 2 (3(2/1 + z/2) 2 + y 2 )) ■ 

Thus if cr > 1 we have V < 0 and the chain of FHN neurons synchronizes. 


□ 


Generalizing two coupled systems to a directed chain of n oscillators, we observe that the Laplacian 
matrix of this configuration is 


0 

0 

0 

•• 0\ 

1 

1 

0 

•• 0 

0 

-1 

1 





0 

0 


0 

-11/ 


The matrix L has only real eigenvalues; a simple zero eigenvalue and n— 1 eigenvalues equal to 1. The only 
type of stable correlated oscillations we can find in the chain are the completely synchronous oscillations. 
These synchronous oscillations will emerge for values of the coupling strength a for which the chain of 2 
FHN oscillators synchronize. Thus the conditions for synchronization are independent of the size of the 
network (i.e. the length of the chain). Numerical simulations below illustrate this statement. 

Figure [2] shows the outputs of two FHN oscillators and the synchronization output error for a = 1.5. 
Figure [3] shows the results for longer chains; Even though the convergence to the synchronous state is 
slower for longer chains, the oscillators in the chains always end up in synchrony. 


3.2.2. Directed ring: ”looking back” 

Suppose now that the n-th oscillator is feeding back its output to the input of the 1st, that is the network 
topology is that of the directed ring. As we shall see later the presence of such an extra connection has 
a drastic effect on the system’s performance with respect to the coupling strength needed to maintain 
stable full-state synchrony. This is reflected in the statement of the theorem below. 
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(a) n = 10 
41-.-i- 



■4 -‘- 1 - 1 - 1 - 

0 200 400 600 800 1000 

t 

(c) n = 100 



(b) n = 50 

4 |-■--- 



-2 


4'- ‘ -'- 1 - 1 - 

0 200 400 600 800 1000 

t 

(d) n = 150 


Figure 3. Synchronization output errors yj := yj — yj- 2 , j = 2, ..., n, for a = 1.5 and 
different lengths of the chain. 


Theorem 3.4. Consider the system of coupled FHN oscillators /TO ). fOD . and let Xi < X 2 <■■■< X n 
be the eigenvalues of the symmetrized Laplacian of the network ^(L + L T ). Then solutions of the coupled 
system asymptotically synchronize providing that 

crX 2 > 1 

Proof. Consider the new variables 

z = Lz, y = Ly, 

where L the Laplacian matrix of the ring, i.e. 


(f) 

Z2 


/ Z\- z n \ 

-2 - 


(y}\ 

V2 


/ Vl-Vn \ 
2/2 - 2/1 

\z n I 


\^n-l - Z n ) 

and y = 

\Vn / 


\Vn-l ~ 2 In) 


It is clear that the systems are synchronized if and only if 

z = 0 and y = 0. 
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Observe that 1 ^ range(-L), hence there exist no vectors z and y such that 

Lz = 1 and Ly = 1. 

This means that the projections of (z,y) via L take values in the set 

n:={(~z,y)eR 2n \z±l,y±l}. 

Thus all synchronization errors 5 and y are orthogonal to 1. 

Consider the function V : £2 —> R + : 

V=\^z T ~z + y T y). 


From the discussion on synchronization in the chain it follows that 

V < ~/3z T z + y T (I — aL)y — y T Wy 


where 


W = 


/ (3(z/i + Vn) 2 + yl) 


\ 


(3(j/n+l + Vn) 2 + Hi), 


which is positive semi-definite, hence 


V < —/3z t z + y T (I — aL)y 

For all vectors yl 1 the following inequality holds true: 

y T (aL - I)y = y T (a\(L + L T ) - I)y > (aX 2 - 1 )y T y 

where A 2 = A 2 (\(L + L T )) is the smallest non-zero eigenvalue of -/(L + L T ). An application of LaSalle’s 
invariance principle, cf. |22|, implies that the synchronization errors z and y converge to zero asymptot¬ 
ically. □ 

Corollary 3.5. For the network of n coupled FHN oscillators, solutions globally asymptotically synchro¬ 
nize if the following inequality holds: 


a 


1 — cos 



> 1 


Proof. Note that 1(L + L T ) is the Laplacian matrix of the undirected ring, which has a simple zero 
eigenvalue with corresponding eigenvector in span(l). According to the properties of the spectrum of 
circulant matrices, cf. [8|, we know that the second smallest eigenvalue A 2 of the symmetrized Laplacian 
i(L + L t ) equals the real part of the smallest (in absolute value) non-zero eigenvalue of L , which we 
denote as 9?(A2(L)). Then if 

+ L t )) = cr3fJ(A 2 (L)) > 1, 

we have V < 0, i.e. V is a Lyapunov function on 17. Note that 


A 2 (A) 


1 


2ir«(r»-l) 

— e 


from which the result immediately follows. 


□ 
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3.3. Synchronization and rotating waves 

The results in the previous sections show that, on the one hand, when a system has a directed ring 
topology and the number of systems in the ring grows then their relative dynamics becomes more and 
more underdamped (Theorem 12.11) . On the other hand, in accordance with Corollary 13.51 estimates of 
attraction rates of the diagonal synchronization manifold rapidly diminish to zero with increasing numbers 
of systems. The latter result is, however, sufficient and may be conservative. To get a clearer view of the 
network dynamics we performed an exhaustive numerical exploration of the system dynamics for various 
values of coupling strengths <r as well as network sizes n. 

We construct a grid (n, a) for number of systems n = 2,..., 20 and coupling strengths a = 
{0.05, 0.1, 0.15,..., 10}. For each (n,<r) 100 sets of initial conditions are drawn uniformly randomly 
from the domain |y#0)| < |-\/3 and |z,;(0)| < /y-\/3, which can be shown to be positively invariant for 
both connectivity configurations (i.e. the directed simple cycle and the directed chain). The MATLAB 
numerical solver ode /5 was used with relative and absolute error tolerances of order 1CD 5 to integrate 
dynamics for a maximum of 20,000 time steps. At regular intervals of 1000, 2000,..., 20, 000 time steps 
we interrupt integration to check for synchronization or rotating wave solutions. After 20,000 time steps, 
if neither synchronization nor a rotating wave solution is detected, we register ‘no solution’. 

Synchronization is identified in terms of the absolute error between the states of neighbouring systems 
averaged over a 1000 time step window being less than 2 x 10 -5 . In case of no synchronization, we 
investigate the existence of rotating waves of Mode Type 1. Rotating waves are defined as periodic 
solutions where all systems take identical orbits with constant non-zero and equal phase shifts between 
neighbouring systems. The mode type describes the group velocity of the wave; for a periodic wave, Mode 
Type 1 describes the case where the period of a rotating wave having non-zero wave velocity equals the 
period of individual oscillators. Identical orbits are identified if the absolute difference between the time 
shifted orbits - so that orbits are in-phase - of neighbouring systems averaged over the period of the orbit 
is less than 10 -4 . Constant and equal phase shifts (for a Mode Type 1 rotating wave) are identified if the 
maximum from all absolute differences between n times the phase shifts between pairwise neighbouring 
systems and period T is less than a tolerance of 10~ 2 . 

The results of this exploration are summarized in Figured This figure shows that in addition to regions 
corresponding to mere full asymptotic synchronization there is a wide range of parameter combinations 
(growing with system size) for which the system admits an asymptotically stable rotating wave solution. 
The larger the number of systems, the larger values of the coupling parameter a are required to maintain 
global stability of the fully synchronous state. 

Two solid curves approximate boundaries between the parameter domains corresponding to analytically 
determined globally asymptotically stable full-state synchrony, and a partition of numerically determined 
globally asymptotically stable full-sate synchrony; The first region corresponds to synchronization reg¬ 
istered for every set of random initial conditions during numerical simulations, whilst the second region 
corresponds to, in addition to synchronization registered for every set of random initial conditions during 
numerical simulations, where Floquet stability analysis of solutions of the auxiliary system indicated 
existence of a locally asymptotically stable rotating wave solution. 

The first (lower) curve - separating analytical and numerical synchronization - was determined previ¬ 
ously in the semi-passivity argument. The second (upper) curve - partitioning numerically determined 
globally asymptomatic synchronization - is determined from a local stability analysis (using Floquet 
theory) of the rotating wave solution. Details of the second are provided below. 

Figure [5] shows for each a and n the proportion of initial conditions that yield a rotating wave solution 
of Mode Type 1 whilst Figure [6] shows for all mode types, i.e. rotating waves that resonate with individual 
systems period of oscillation. For low coupling a and for increasing number of systems n, rotating wave 
solutions are found more often. This suggests a larger basin of attraction for the rotating wave than that 
for synchronization, and that this basin grows with increasing n and decreasing a whilst at the same 
time the basin of attraction for synchronization shrinks. The relative sizes of basin of attraction result 
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Figure 4. Bifurcation diagram for directed rings of FHN oscillators. The diagram is 
divided into the four regions of parameter space corresponding to: Synchronization (1), 
global asymptotic synchronization that is guaranteed by the semi-passivity argument; 
Synchronization (2), synchronization registered for every set of random initial condi¬ 
tions during numerical simulations; Synchronization (3) synchronization registered for 
every set of random initial conditions during numerical simulations, but Floquet stability 
analysis of solutions of the auxiliary system indicated existence of a locally asymptoti¬ 
cally stable rotating wave solution; Co-existence , both the fully synchronous and rotating 
wave solutions were registered during numerical simulation. 


in higher or lower likelihoods for the systems to converge to a certain solution given uniformly random 
initial conditions. 


3.3.1. Local stability analysis of the rotating wave 


Throughout this section we consider only the rotating wave of Mode Type 1. Similar analysis can also 
be performed for other mode types. 


Suppose that n identical coupled systems have a non-constant, T-periodic 
constant T > 0, and for which the orbit of each system is identical and time 



solution Xj = ( Zj,yj ) for 
shifted by some constant 


xi (t) = x 2 (t + t) = x 3 (t + 2t) = ■ ■ ■ = x n (t + (n — 1 )t) 

= Xi(t + nr) = X\(t + T). (3.7) 

We refer to this as the rotating wave solution. An example of a rotating wave solution for n = 5 coupled 
FHN oscillators in the ring configuration is presented in Figure [7] 

Recall equation (13.11) with Xj = ( Zj,yj). If we restrict the coupled dynamics of the FHN oscillators 
to the rotating wave manifold, then using the periodicity of the rotating wave solution, substitution of 
equation (13.71) into the dynamics of each coupled FHN oscillator (13.11) yields n identical uncoupled delay 
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Figure 5. Proportion of samples that yield a rotating wave solution of Mode Type 1. 
Red curve bounds the upper left region for which Floquet stability analysis indicated the 
existence of a rotating wave solution. 



0.05 12 34 56789 
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Figure 6. Proportion of samples that yield a rotating wave solution for all mode types. 


differential equations (DDEs) of the form 

xi (i) = f(x i (i)) + crBC (t - - x\ (t)j 


x n (t) = f(x n (t)) + aBC 



^^T)-X n (t)). 
n J 
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T 




(a) (b) 

Figure 7. (a) The y dynamics over one oscillation, (b) The periodic y dynamics in the 
time interval [0, T]. The z dynamics show the same type of time shifted and periodic 
behavior as the y dynamics. 


Thus the rotating wave solution can only exist if the auxiliary system 

s(t) =/(s(0) - aBC[s(t) - s(t - r*)], 

* rp n ~ l rr 

t :=1 — t = -7, 

n 

has a non-constant, T-periodic solution: 

St = St+T £ C = C([0, T], K 2 ), 

for which the set C is the set of continuous functions that map the interval [0,7/ into 

s(t + d ), 9 £ [0, T\. 

Define the errors between neighbouring systems around the rotating wave solution 

ej(t) = x j+1 (t + r) - Xj(t), j = l,2,...,n- 1, 

e n (t) = X!(t + t) - x n (t). 

Taking the error dynamics we obtain 


(3.8) 

(3.9) 
and s t {0) := 

(3.10) 


/ei(0\ 

6(0 


( /(ei(0 +Xi(0) - fi x l(O) \ 

/(e 2 (0 +x 2 {t)) - f {x 2 (0) 

- cr{L ® BC) 

( e i(0 \ 

e 2 (0 

\ e-nit) ) 


V/(e n (0 + x n{t)) - f{x n {t))j 


\ e n {t) J 


(3.11) 


Substitution of the rotating wave solution in terms of the auxiliary system variable s{t) into equation 
(13.111) . such that 

s(t) = x\(t) = xz(t + t) = • • • = x n (t + (n - l)r), 
and linearizing around the rotating wave solution yields the linear system (13.121) 


(i i(*)\ 

6(0 

\Ut)J 


/ J(s(0) 




[\ 


\ 

J(s(t - (n- 1 )t)) / 
15 


— a(L (8) BC) 


/ 6 ( 0 \ 

6(0 


(3.12) 


. \ 6 ( 0 / 
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Figure 8. Solutions of the auxiliary system characterized in the parameter domain of 
period time, delay, and coupling strength (T, r, a) presented as a surface for T a function 
of pairs (r, a). 


where ® is the Kronecker (tensor) product, J(s(t )) is defined as follows: 


W»:= ('“/ !-“*))■ 

and S 2 {t) denotes the second component of s(t). Note that T-periodicity of the system (13.81) implies the 
linear error (13.121) system to be T-periodic. 

For the local stability analysis we first computed periodic solutions of the auxiliary system (13.81) . 
Periodic solutions are determined using continuation methods that are available in the numerical software 
package DDE-Biftool jlj |. Figure (JSJ) characterizes solutions of the auxiliary system in the parameter 
domain (T, r, a). For the auxiliary system in which parameter T and r are allowed to vary continuously, 
a solution that describes the dynamics of a rotating wave solution satisfies the relation -^(n — 1) = n. 
However, for the solutions we obtained, parameters T and r have not been varied continuously. Therefore, 
we choose the solution that satisfies the following the inequality 


T 

— (n — 1) — n 

T 


< e. 


(3.13) 


To maintain good accuracy of approximation of the auxiliary system to n coupled FHN oscillators, the 
error e must be small. For our stability analysis we took e = 0.01. 

Figures 9(a) and |9(b)| show two cross sections of the surface in Figurc[8]for coupling strengths a = 0.95 
and a = 6.75, respectively. Dashed lines identify solutions that satisfy relation (13.131) and hence map 
solutions of the auxiliary system to an integer number n of coupled FHN oscillators on the rotating wave 
manifold. 

We assessed the stability of the rotating wave solution for pairs (n, a) by computing the Floquet 
multipliers of the periodic linearized error system (13.121) (again with DDE-Biftool). Solutions are obtained 
by substituting in the solution of the auxiliary system corresponding to the pair (n, a) found by numerical 
continuation. Recall that if all Floquet multipliers except one (at 1) have modulus strictly smaller than 1, 
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Figure 9. Solutions of the auxiliary system in the parameter domain of period time 
and delay (T,r) for given coupling strengths: (a) a = 0.95. (b) cr = 6.75. Dashed lines 
indicate solutions that satisfy relation (13.131) . 


then the zero solution of the linearized error system is asymptotically stable, which implies the rotating 
wave solution to be locally orbitally stable. The red line in Figure []] (and Figure [5]) is defined by the 
crossing of (at least) one multiplier with the boundary of the unit disc in C. 


4. Discussion 

4.1. Kinetic interpretation of (12.11) 

Equation m in Section [2] describes the temporal evolution of the first order kinetics. This equation 
is known as the Master Equation. The master equation obeys the principle of detailed balance if there 
exists a positive equilibrium P* (jp* > 0) such that for each pair i,j (i # j) 


qijPj = q 0 iPi ■ 


(4.1) 


After Onsager [24], it is well known that for systems with detailed balance the eigenvalues of K are real 
because under conditions (EH) IF is a self-adjoined matrix with respect to the entropic inner product 


(x, y)=Yl 


Xiyi 

P* 


(see, for example, |3lll34l ]). 

Detailed balance is a well known consequence of microreversibility. This principle was introduced in 
1872 by Boltzmann for collisions [§!]. In 1901 Wegscheider proposed it for chemical kinetics j33|. Einstein 
had used it as a principle for the quantum theory of light emission and absorption (1916, 1917). The 
backgrounds of detailed balance had been analyzed by Tolman [30j ]. The principle was studied further 
and generalized by several authors mum. 


Systems without detailed balance appear in applications rather often. Usually, they represent a sub¬ 
system of a larger system, where concentrations of some of the components are considered as constant. 
For example, the simple cycle 

A\ —> A 2 A n —> A\ (4-2) 
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is a typical subsystem of a catalytic reaction (a catalytic cycle). The complete reaction may have the 
form 

S + A\ —> A 2 —y ... —y A n —y A\ + P , (4-3) 

where S' is a substrate and P is a product of reaction. 

The irreversible cycle (14.21) cannot appear as a limit of systems with detailed balance when some of 
the constants tend to zero, whereas the whole catalytic reaction (14.31) can [17]]. The simple cycle (14.21) 
can be produced from the whole reaction (14.31) if we assume that concentrations of S and P are constant. 
This is possible in an open system, where we continually add the substrate and remove the product. 
Another situation when such an approximation makes sense is a significant excess of substrate in the 
system, [5] [Ai\ (here we use the square brackets for the amount of the component in the system). 
Such excess implies separation of time and the system of intermediates {A^} relaxes much faster than 
the concentration of substrate changes. 

In systems without detailed balance, damped oscillations are possible. The example in Section [2j which 
describes the case of all the reaction rate constants in the simple cycle being equal, Qj+ij = q\ n = q > 0, 
shows that these oscillations are even weakly damped. The effect becomes acutely noticeable for n large 
enough. 

The simple cycle with equal rate constants yields the slowest decay of oscillations or, in some sense, 
the slowest relaxation among all first order kinetic systems with the same number of components. The 
extremal properties of the simple cycle with equal constants were noticed in numerical experiments 25 
years ago Q. V.I. Bykov formulated the hypothesis that this system has extremal spectral properties. 
This paper provides the answer: yes, it has. 


4.2. Two coupled cycles 

Given the size of the region where multiple solutions co-exist, and the resilience to a coherent state; does 
the extremal property of the simple cycle give rise to further, more complex phenomena when two simple 
cycles are diffusively coupled via an undirected link between an oscillator in each cycle? 

For a total of 2k coupled systems, two cycles are constructed with systems 1,..., k in the first simple 
cycle and systems k + 1,... ,2k in the second, and coupled via systems X\ and Xk+i- Clearly the syn¬ 
chronization manifold exists, as does the rotating wave solution in the form of two synchronized rotating 
waves, 


xi (t) =x k +i{t) = x 2 {t + r) = x k+2 (t + r) = ... 

■ ■ ■ = x k -i(t + (k- 2)t) = x 2 k-i(t + (k - 2)t ) = Xk(t + [k - l)r) = x 2k (t + (k - l)r). 

A full description of the phenomena of two coupled cycles is beyond the scope of this work; however, as 
a motivation for further study, we present a brief example. 

We take (n, a) = (10, 0.75), which, for a simple cycle lies in the region of co-existence of synchronization 
and rotating wave solutions. We observe in Figure (fTOl) a stable state in which the trajectories of all 
systems in the first cycle (in red) are attracted to the synchronization manifold, whilst all trajectories 
of systems in the second cycle (in green) are attracted to the rotating wave solution. There is a clear 
competition of each cycle to attract the other to its own dynamical regime. The two diffusively coupled 
oscillators from each cycle periodically perturb each other, which prevents asymptotic convergence of 
systems to either the synchronization manifold or the rotating wave solution. Clearly, the extremal 
properties of the simple cycle can give rise to multiple regimes of complex patterns of dynamics when 
embedded into larger network structures. 

5. Conclusion 

We considered the problem of how “closing” a chain of interconnected systems with directed coupling by 
adding a directed feedback from the last element in the chain to the first may affect the dynamics of the 
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Figure 10. Two coupled cycles and their y-dynamics; in red the ^-dynamics of the 
first cycle, and in green the ^-dynamics of the second cycle. 


system. This problem is closely related to the fundamental question of how network topology influences 
the dynamics of collective behavior in the system. Two general settings have been investigated. In the 
first one we analyzed the behavior of a simple linear system. We showed that the simple cycle with 
equal interaction weights has the slowest decay of the oscillations among all linear systems with the same 
number of states. In the second setting we considered directed rings and chains of identical nonlinear 
oscillators. For directed rings, a lower bound a c for the connection strengths that guarantee asymptotic 
synchronization in the network is found to follow a pattern similar to that of a simple cycle. Furthermore, 
numerical analysis revealed that, depending on the network size n, multiple dynamic regimes co-exist in 
the system’s state space. 


In addition to the fully synchronous state, for sufficiently large networks an asymptotically stable 
rotating wave solution emerges. The emergence of the rotating wave is a phenomenon that persists 
over a broad range of coupling strengths and network sizes, and can be viewed as a form of extreme 
sensitivity of the network dynamics to the removal or addition of a single connection. The result confirms 
the significance of shortcuts in networks with large numbers of nodes. Emergence of asymptotically 
stable rotative wave solutions has been analyzed numerically for a specific class of systems in which 
the dynamics of each node was identical and satisfied Fitzhugh-Nagumo equations 11]. Extending the 
analysis to systems with heterogeneous nodes as well as considering nodes with Hindmarsh-Rose and 
Hodgkin-Huxley dynamics []j|, known to be capable of bursting and chaotic behavior, will be the topic 
of our future studies. 


Coming back to the question if leaders should look back. To stay in synchrony we advise a leader 
either not to look back at all or to look back just a few links; looking back too far induces oscillations 
that destroy the coherent state. 
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